library(magrittr)
x <- 9
# Calculate the square-root of x
sqrt(x)

## [1] 3

# Calculate it using pipes
x %>% sqrt

## [1] 3

x <- 9
# Calculate the square-root of x and update x
x <- sqrt(x)
x

## [1] 3

# Calculate it using pipes
x <- 9
x %<>% sqrt
x

## [1] 3

df <- read.csv("beedata.csv")

nrow(subset(df, hive==4))

## [1] 60193

df %>% subset(hive==4) %>% nrow

## [1] 60193

# Exercises - Piping

# 1. Rewrite the following code using %>% and %<>%:

x <- 2
round(log(x))

## [1] 1

# 2. Rewrite the second line of following code:

x <- rnorm(10,100)
round(sum(sqrt(x)), 3)

## [1] 100.027

x <- 2
round(log(x))

## [1] 1

x <- 2
x %>% log %>% round

## [1] 1

# 2. Rewrite the second line of following code:

x <- rnorm(10,100)
round(sum(sqrt(x)), 3)

## [1] 99.993

x %>% sqrt %>% sum %>% round(3)

## [1] 99.993

Artwork by @allison_horst

Horst AM, Hill AP, Gorman KB (2020). palmerpenguins: Palmer Archipelago (Antarctica) penguin data. R package version 0.1.0. https://allisonhorst.github.io/palmerpenguins/. doi: 10.5281/zenodo.3960218.

https://allisonhorst.github.io/palmerpenguins/

library(palmerpenguins)

## 
## Attaching package: 'palmerpenguins'

## The following objects are masked from 'package:datasets':
## 
##     penguins, penguins_raw

head(penguins)

## # A tibble: 6 × 8
##   species island    bill_length_mm bill_depth_mm flipper_length_mm body_mass_g
##   <fct>   <fct>              <dbl>         <dbl>             <int>       <int>
## 1 Adelie  Torgersen           39.1          18.7               181        3750
## 2 Adelie  Torgersen           39.5          17.4               186        3800
## 3 Adelie  Torgersen           40.3          18                 195        3250
## 4 Adelie  Torgersen           NA            NA                  NA          NA
## 5 Adelie  Torgersen           36.7          19.3               193        3450
## 6 Adelie  Torgersen           39.3          20.6               190        3650
## # ℹ 2 more variables: sex <fct>, year <int>

https://github.com/mcnakhaee/palmerpenguins?tab=readme-ov-file

from palmerpenguins import load_penguins

## /home/diseng001/R/x86_64-pc-linux-gnu-library/4.5/reticulate/python/rpytools/loader.py:120: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.
##   return _find_and_load(name, import_)

penguins = load_penguins()
penguins.head()

##   species     island  bill_length_mm  ...  body_mass_g     sex  year
## 0  Adelie  Torgersen            39.1  ...       3750.0    male  2007
## 1  Adelie  Torgersen            39.5  ...       3800.0  female  2007
## 2  Adelie  Torgersen            40.3  ...       3250.0  female  2007
## 3  Adelie  Torgersen             NaN  ...          NaN     NaN  2007
## 4  Adelie  Torgersen            36.7  ...       3450.0  female  2007
## 
## [5 rows x 8 columns]

• What is your response variable?
• Plot your data
  • Plot raw data
  • Frequency distribution -> Are there any patterns?

\[ \begin{equation} y_i = \beta_0 + \beta_1x_i + \epsilon_i, \qquad i=1,...,n \end{equation} \] where \(y\) is the dependent variable, \(x\) is the independent variable (also called explanatory variable), \(\beta_0\) is the intercept parameter, \(\beta_1\) is the slope parameter and \(\epsilon\sim N(0,\sigma^2)\) is the error coefficient.

set.seed(123)
x <- seq(0,5,0.1) 
y <- x + rnorm(length(x))
plot(x, y)

mod <- lm(y~x)
summary(mod)

## 
## Call:
## lm(formula = y ~ x)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.0116 -0.6110 -0.0912  0.6575  2.1444 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.05829    0.25565   0.228    0.821    
## x            0.99216    0.08812  11.259  3.4e-15 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9263 on 49 degrees of freedom
## Multiple R-squared:  0.7212, Adjusted R-squared:  0.7155 
## F-statistic: 126.8 on 1 and 49 DF,  p-value: 3.397e-15

plot(x,y)
abline(reg = mod)

Our hypothesis is that there is no relationship between the variables. The p-value expresses the probability of this hypothesis. If our p-value is very small, it is very unlikely, that this hypothesis is true and we have a statistical significant relationship between the variables.

• Linear relationship between x and y
• Independence of residuals
• Constant variance of residuals
• Normality: Residuals are normally distributed

plot(mod, which=1)

• 1st graph: Constancy of variance
• Plot of the residuals against the fitted values
• Should look like the sky at night

plot(mod, which=2)

• 2nd graph: Normality of errors

• Normal QQ plot = plot of the ordered residuals against the fitted values

• Should be a reasonably straight line

plot(mod, which=3)

plot(mod, which=4)

• use a scatterplot to check
• if this assumption is violated we have unrealistic predictions and estimates of the standard deviation
• How to deal with this? You can transform the variables or alter the model.

plot(x,y)

• Plot residuals in time/ observation order / based on observation location / …
• You expect that they are equally scattered around 0
• You can also use the Durban-Watson-Test
• Autocorrelation can have serious effects on the model
• You can add an AR correlation term to the model. Sometimes it also helps to add a further covariate.

res <- residuals(mod)
plot(x, res)

• H0: There is no correlation among the residuals.

• HA: The residuals are autocorrelated.

library(car)

## Loading required package: carData

durbinWatsonTest(mod)

##  lag Autocorrelation D-W Statistic p-value
##    1      0.02431999      1.940949   0.712
##  Alternative hypothesis: rho != 0

• The p-value is larger than 0.05, so we cannot reject the null hypothesis that there is no correlation between residuals.

• Plot fitted values vs residuals
• You expect that they are equally scattered around 0

plot(fitted(mod),res)

• Use qq-plot to compare quantiles.
• Or use the Shapiro-Wilk-Test
• If you do not find a Normal Distribution, check for outliers or transform your variable.

qqPlot(res)

## [1] 44 18

• H0: The data is normally distributed
• HA: The data comes from an other distribution

shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.99152, p-value = 0.9732

• We cannot reject the null-hypothesis that the data is normally distributed.

If not another model might be more appropriate:

• Continuous: linear model → LM e.g. height, weight
• Count: Poisson model → GLM e.g. number of individuals
• Binary: Binomial model → GLM e.g. presence / absence

For example many zeros:

You could separate in two analyses: Binomial model for success / failure and linear model for successes

4th graph of plot(model): Plot of Cook´s distances versus raw data → Highlights the identity of particularly influential data points

Polynomial regression e.g. a quadratic function might be an option

 lm(y ~ x + I(x^2))

You could for example use a log transformation of the response variable

    lm(log(y) ~ x)

\[\begin{equation}\label{eqn:linearregression} y_i = \beta_0 + \beta_1x_{1,i} + ... + \beta_px_{p,i} + \epsilon_i, \qquad i=1,...,n \end{equation}\] where \(y\) is the dependent variable, \(x_1 ... x_p\) are the independent variables (also called explanatory variables), \(\beta_0 ... \beta_p\) are the regression coefficients, \(\epsilon\sim N(0,\sigma^2)\) is the error coefficient and \(p \geq 1\).

• many, many libraries
• we are going to look at base, ggplot and plotly

df.4 <- df[df$hive==4,]
# plot temperature outside
plot(df.4$time, df.4$t_i_3)

# plot temperature outside
plot(df.4$time, df.4$t_i_3, 
     ylim=c(0,40)) # added this

# plot temperature outside
plot(df.4$time, df.4$t_i_3, 
     ylim = c(0,40), 
     xlab = "Time (2019)",  # added this
     ylab = "Temperature within hive", # added this
     main = "Sensor measurements") # added this

# plot temperature outside
plot(df.4$time, df.4$t_i_3, 
     ylim = c(0,40), 
     type = "b", # added this
     lty = 1, # added this
     xlab = "Time (2019)", 
     ylab = "Temperature within hive", 
     main = "Sensor measurements")

• “p”: Points
• “l”: Lines
• “b”: Both

Source: http://www.sthda.com/english/wiki/line-types-in-r-lty

df.4 <- df.4[df.4$t_i_3>5&df.4$t_i_3<40,] # added this
# plot temperature outside
plot(df.4$time, df.4$t_i_3, 
     ylim = c(0,40), 
     type = "b",
     lty = 1,
     xlab = "Time (2019)", 
     ylab = "Temperature within hive", 
     main = "Sensor measurements")

# plot temperature outside
plot(df.4$time, df.4$t_i_3, 
     ylim = c(0,40), 
     type = "b",
     lty = 1,
     pch = 4, # added this
     xlab = "Time (2019)", 
     ylab = "Temperature within hive", 
     main = "Sensor measurements")

Source: http://www.sthda.com/english/wiki/r-plot-pch-symbols-the-different-point-shapes-available-in-r

# plot temperature outside
plot(df.4$time, df.4$t_i_3, 
     ylim = c(0,40), 
     type = "b",
     lty = 1,
     pch = 4,
     xlim = as.POSIXct(c("2019-08-08", "2019-08-09")), # added this
     xlab = "Time (2019-08-08)", 
     ylab = "Temperature within hive", 
     main = "Sensor measurements")

# plot temperature outside
plot(df.4$time, df.4$t_i_3, 
     ylim = c(0,40), 
     type = "b",
     lty = 1,
     pch = 4,
     xlim = as.POSIXct(c("2019-08-08", "2019-08-09")),
     xlab = "Time (2019-08-08)", 
     ylab = "Temperature within hive", 
     main = "Sensor measurements",
     xaxt="n")
axis.POSIXct(1, 
             at=seq(min(df.4$time), max(df.4$time), by="1 hour"), 
             format="%H:00") # added this

# subset data
df.4 <- df[df$hive==4,]
# plot temperature outside
plot(df.4$time, df.4$t_o, ylim=c(0,40),type = 'p', pch=4)
# choose colours
cl <- rainbow(5)
# choose colums
cols <- 4:8
# plot each column
for (i in 1:5){
    lines(df.4$time, 
          df.4[,cols[i]],
          col = cl[i],
          type = 'p', 
          pch=4, 
          ylim=c(0,40))
}
# add legend
legend("topright", legend=c(1, 2, 3, 4, 5, "outside"),
       col=c(cl, "black"), pch = 4, lty = 0, cex=0.8)

# add legend
legend("topright", legend=c(1, 2, 3, 4, 5, "outside"),
       col=c(cl, "black"), pch = 4, lty = 0, cex=0.8)

# plot data
library(ggplot2)
ggplot(data = df.4, aes(x=time, y=t_i_3)) + geom_point()

# plot data
library(ggplot2)
ggplot(data = df.4, aes(x=time, y=t_i_3)) + geom_point(shape=4) + 
  ylim(c(0, 40)) + 
  xlab("Time (2019") + 
  ylab("Temperature within hive") + 
  ggtitle("Sensor measurements")

## Warning: Removed 152 rows containing missing values or values outside the scale range
## (`geom_point()`).

# subset data
df.4 <- df[df$hive==4,]
# choose columns
df.4.cols <- df.4[,c(1,4:9)]
# reshape data
library(reshape)
mdf <- melt(df.4.cols, id=c("time")) 
# plot data
library(ggplot2)
ggplot(data = mdf, aes(x=time, y=value)) + 
  geom_line(aes(colour=variable)) + 
  ylim(c(0, 40))

library(plotly)
fig <- plot_ly(df.4[1:100,], x = ~time, y = ~t_i_3)
fig

p <- plot_ly(df.4[1:100,], x = ~time, y = ~t_i_3)
# saving the plot as html
htmlwidgets::saveWidget(p, "testploty.html")

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plotly.com/r/reference/#scatter

## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode

# if you want to have a static image, e.g. png, you can then open the html file, and use the small button with the camera symbol to export as png. You could also use the library orca to directly export a static image using plotly directly from your code, but it seems to be a bit complicated to install orca.
# I haven't used powerpoint for a long time, but I think there should be a way to embed a html file: https://www.techwalla.com/articles/how-to-embed-an-excel-workbook-icon-into-powerpoint

png("test.png")
plot(hist(rnorm(100)))
dev.off()

• you all get a confirmation of participation
• you can submit a small assignment for 1 CP

• you can choose if you want to use R or Python
• Create a markdown document/ Jupyther Notebook
• Choose a dataset from PANGEA
• Think about the following question: What do you want fo find out? What is the motivation for your analysis?
• Choose 5 of the following tasks. You can also come up with own ideas.
• Write a short text for each task explaining what you have been doing

• Make some exploratory analysis:
  • print the mean/ median/ standard deviation for each column
  • print the total number of missing values
  • print the number of missing values for each column
• Transform some of the columns/ Create new columns based on existing ones
• Create a subset of your data
• Exclude all rows with missing values
• Create a plot
• Write a function that takes a filename and some additional parameters. The function should create a plot and save it as “png” with that filename
• Create a linear model

• If you want to do an assignment, please send me an e-mail:

• “I want to do an assignment.””

• “I do not want a grade on my certificate. / I want to have a grade on my certificate./ I want to have a grade if it is better than <1.3, 1.7, 2.0, 2.3, 2.7, 3.0, …>”

• Suggestion for deadline: ?